This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/ERL_Angourakis-et-al-2020_Rproject

Figure 3. Daily precipitation and seasonal totals (mm) for on the location of the five major Indus urban centres

Using:

Daniel Wallach, David Makowski, James W. Jones, François Brun (2019), ‘Working with dynamic crop models: Methods, tools, and examples for agriculture and enviromnent’, Academic Press.

Model description in p. 24-28, R code example in p. 116-122. Another, possibly newer version of the R code can be found at: https://github.com/cran/ZeBook/blob/master/R/watbal.model.r

Original model from:

Woli P, Jones J W, Ingram K T and Fraisse C W (2012) ‘Agricultural Reference Index for Drought (ARID)’ Agron. J. 104-287. Online: https://www.agronomy.org/publications/aj/abstracts/104/2/287

Load functions from ‘Soil Water Balance model’ library file:

source("source/watbal.model.R")

Technical note: The ZeBook package, watbal.model.R has a bug in lines 21-25 and 58. The problem is that the code fails to extract values from the param object generated by watbal.define.param because it does not account for both dimensions of this object: (WHC, MUF, DC, z, CN) x (nominal, binf, bsup). I solve this problem by adding “typeOfParameterValue” as an argument in watbal.model that is passed also to watbal.update.

Load function to estimate reference evapotranspiration:

source("source/estimateETr.R")

Load CurlyBraces function for plotting:

source("source/CurlyBraces.R")

< begin parentesis >

To reproduce Figure 4.5 in Wallach et al. 2006, load weather daily data and select site and year :

weather <- watbal.weather(working.year = 2008, working.site = 1)

Or following code on the book:

weather <- ZeBook::weather_FranceWest
weather <- weather[(weather$idsite == 1 & weather$WEYR == 2008),]

And skip to the steps related to the Soil Water Balance model (bellow).

< end parentesis >

To generate Figure 3, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ and pin pointing the five different locations between 01/01/1984 and 31/12/2018. The exact locations are:

We selected the ICASA Format’s parameters:

The following code reads all files inside the “Fig3_input” folder and creates a single data frame that contains the data of all five locations:

inputFiles <- paste0("Fig3_input/", list.files(path = "Fig3_input"))

weather <- data.frame()

for (i in 1:length(inputFiles))
{
  tempData <- watbal.weather.file(inputFiles[i], year = 1984:2008)
  
  # convert any missing data in all sky insolation (I) from -99.0 to Na
  tempData$I[tempData$I == -99] <- NA
  
  weather <- rbind(weather, tempData)
}

rm(tempData)

# use Latitude to distinguish sites and create a new variable using site names:
weather$Site <- rep(NA, nrow(weather))
for (i in 1:nrow(weather))
{
  if (floor(weather$LAT[i]) == 27) { weather$Site[i] <- "Mohenjo-daro" }
  if (floor(weather$LAT[i]) == 23) { weather$Site[i] <- "Dholavira" }
  if (floor(weather$LAT[i]) == 28) { weather$Site[i] <- "Ganweriwala" }
  if (floor(weather$LAT[i]) == 30) { weather$Site[i] <- "Harappa" }
  if (floor(weather$LAT[i]) == 29) { weather$Site[i] <- "Rakhigarhi" }
}
weather$Site <- factor(weather$Site, 
                       levels = c("Mohenjo-daro", "Ganweriwala", "Harappa",
                                  "Dholavira", "Rakhigarhi"))

Sum up yearly and seasonal totals by splitting the year into summer (May-Oct or from day ) and winter (Nov-Apr):

precipitationSums <- list(Site = c(), year = c(), total = c(), summer = c(), winter = c())

for (site in levels(weather$Site))
{
  for (year in levels(factor(weather$year)))
  {
    precipitationSums$Site <- c(precipitationSums$Site, site)
    precipitationSums$year <- c(precipitationSums$year, year)
    
    precipitationSums$total <- c(precipitationSums$total,
                                 round(sum(weather$RAIN[weather$Site == site & 
                                                    weather$year == year])))
    precipitationSums$winter <- c(precipitationSums$winter,
                                  round(sum(weather$RAIN[weather$Site == site & 
                                                     weather$year == year &
                                                     (weather$day <= 121 | weather$day > 305)]
                                      )))
    precipitationSums$summer <- c(precipitationSums$summer,
                                  round(sum(weather$RAIN[weather$Site == site & 
                                                     weather$year == year &
                                                     (weather$day > 121 & weather$day <= 305)]
                                      )))
  }
}

precipitationSums <- data.frame(precipitationSums)

Calculate the average annual precipitation in every site to reference in the text:

for (site in levels(precipitationSums$Site))
{
  cat(paste0(site, "\n"))
  print(summary(precipitationSums$total[precipitationSums$Site == site], na.rm = TRUE))
  cat("---------------------------------------\n")
}
## Dholavira
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    55.0   188.0   258.0   303.3   420.0   671.0 
## ---------------------------------------
## Ganweriwala
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      30      85     107     117     149     219 
## ---------------------------------------
## Harappa
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     123     164     232     240     299     459 
## ---------------------------------------
## Mohenjo-daro
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    8.00   33.00   72.00   78.96   98.00  286.00 
## ---------------------------------------
## Rakhigarhi
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     248     392     495     499     615     740 
## ---------------------------------------

Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):

weather$ETr <- estimateETr(weather$I, 
                           weather$T2M, 
                           weather$Tmax, 
                           weather$Tmin, 
                           weather$T2MDEW, 
                           weather$WS2M)

Use Soil Water Balance model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.

The first step is to set all soil parameters, here assumed to be constants throughout the sites:

parameters <- watbal.define.param()
# This function creates the following configuration (from Wallach et al. 2006):
# WHC = 0.15
# MUF = 0.096
# DC = 0.55
# z = 40
# CN = 65

# input variables describing soil (from Wallach et al. 2006)
soil.WP = 0.06
soil.FC = soil.WP + parameters["nominal", "WHC"]

Run the model inputting parameters and dataset:

aWatbal <- cbind(Site = weather$Site,
                 watbal.model(parameters, 
                        weather, 
                        WP = soil.WP, 
                        FC = soil.FC, 
                        typeOfParValue = "nominal"))

watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.

Plot precipitation (RAIN), evapotranspiration (ETr), soil water content proportion (WATp), and the ARID coefficient for the five sites during the year 1990.

selectedYear = 1990
lastDayOfWinter = 121
lastDayOfSummer = 305

#--------
# To create a png file:
plotName = "Fig3.png"

fontScale = 1.2

png(plotName, width = 1000, height = 600)
#---------
# alternatively, to create eps file:
# plotName = "Fig3.eps"
# 
# fontScale = 1
# 
# extrafont::loadfonts(device = "postscript")
# grDevices::postscript(file = plotName,
#                       pointsize = 10,
#                       width = 10,
#                       height = 6,
#                       horizontal = FALSE,
#                       paper = "special",
#                       onefile = FALSE,
#                       family = "sans",
#                       colormodel = "cmyk")
#---------

layout(matrix(1:25, nrow = 5, ncol = 5), 
       heights = c(0.2, 1.1, 1, 1, 1.3),
       widths = c(1.25, 1, 1, 1, 1))

yLabs <- c("Precipitation (mm)", "ETr (mm)", "WATp (%)", "ARID index")
leftMargin <- 4.1

# x coordinate in barplots has different scale, so must be adjusted with this offset 
curlyBracesOffset = 1.2 

for (site in levels(factor(weather$Site)))
{
  if (site != levels(factor(weather$Site))[1])
  {
    yLabs <- rep("", 4)
    leftMargin <- 2.1
  }
  tempData <- aWatbal[aWatbal$year == selectedYear & aWatbal$Site == site,]
  
  par(mar = c(0, leftMargin, 0, 0), cex = fontScale)
  
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.7, font = 4, cex = 1.2,
      labels = site)
  text(x = 0.5, y = 0.1, cex = 0.8,
       labels = paste0(precipitationSums$total[precipitationSums$Site == site & 
                                                 precipitationSums$year == selectedYear],
                       " mm (annual total)")
       
       )
  
  par(mar = c(1.1, leftMargin, 2.1, 0.2))

  barplot(tempData$RAIN, 
          ylab = yLabs[1])
  CurlyBraces(x0 = (lastDayOfWinter + 1) * curlyBracesOffset,  
              x1 = lastDayOfSummer * curlyBracesOffset,  
              y0 = 1.1 * max(tempData$RAIN), y1 = 1.1 * max(tempData$RAIN), 
              position = 3, 
              bracesSize = 0.15 * max(tempData$RAIN),
              label = paste0(
                precipitationSums$summer[precipitationSums$Site == site & 
                                                 precipitationSums$year == selectedYear], 
                " mm (summer)"),
              labelSize = 1, labelDist = 1.1)
  CurlyBraces(x0 = 1 * curlyBracesOffset,  
              x1 = lastDayOfWinter * curlyBracesOffset,  
              y0 = -0.1 * max(tempData$RAIN), y1 = -0.1 * max(tempData$RAIN), 
              position = 1, 
              bracesSize = 0.15 * max(tempData$RAIN),
              label = paste0(
                precipitationSums$winter[precipitationSums$Site == site & 
                                                 precipitationSums$year == selectedYear], 
                " mm (winter)"),
              labelSize = 1, labelDist = 1, labelAdj = 0)
  CurlyBraces(x0 = (lastDayOfSummer + 1) * curlyBracesOffset, 
              x1 = 365 * curlyBracesOffset, 
              # x1 coordinate must be higher than 365 to cover the rest of the axis 
              y0 = -0.1 * max(tempData$RAIN), y1 = -0.1 * max(tempData$RAIN), 
              position = 1, 
              bracesSize = 0.15 * max(tempData$RAIN))
  par(mar = c(0.1, leftMargin, 1.1, 0.2))
  
  barplot(tempData$ETr, 
          ylab = yLabs[2])
  
  plot(1:nrow(tempData), tempData$WATp * 100, 
     xlab = "", xaxt = 'n',
     ylab = yLabs[3],
     type = "l", lwd = 2, ylim = c(0, soil.FC * 110))
  abline(h = soil.FC * 100, lty = 2)
  abline(h = soil.WP * 100, lty = 2)
  if (site == levels(factor(weather$Site))[1])
  {
    text(x = 5, y = 17.5, labels = "field capacity", 
         cex = fontScale * 0.8, font = 3, adj = 0)
    text(x = 5, y = 4.5, labels = "wilting point", 
         cex = fontScale * 0.8, font = 3, adj = 0)
  }
  
  par(mar = c(3.5, leftMargin, 1.1, 0.2))
  
  plot(1:nrow(tempData), tempData$ARID, 
      xlab = "day\n", ylab = yLabs[4],
      type = "l", lwd = 2, ylim = c(0, 1))
}
rm(tempData)

dev.off()
## png 
##   2
knitr::include_graphics(plotName)

Figure 6. Example of simulated weather variables

Load output generated with NetLogo implementation of the SIMPLE crop model, which integrates the Weather and Soil Water Balance models:

simExample <- read.csv("Fig6_input/SIMPLE-crop-model experiment-table.csv", skip = 6)

wheatHarvestDay = simExample$"item.0.sugHarvestingDay"[1]
riceHarvestDay = simExample$"item.1.sugHarvestingDay"[1]

simExample_pars <- simExample[1, c(
  "solar_annual.min", "solar_annual.max", "solar_daily.mean.fluctuation",
  "temperature_annual.min.at.2m", "temperature_annual.max.at.2m",
  "temperature_daily.mean.fluctuation",
  "temperature_daily.lower.deviation",  "temperature_daily.upper.deviation",
  "precipitation_yearly.mean", "precipitation_yearly.sd",
  "precipitation_daily.cum_plateau.value_yearly.mean",
  "precipitation_daily.cum_inflection1_yearly.mean", "precipitation_daily.cum_inflection1_yearly.sd",
  "precipitation_daily.cum_rate1_yearly.mean",  "precipitation_daily.cum_rate1_yearly.sd",
  "precipitation_daily.cum_plateau.value_yearly.sd", "precipitation_daily.cum_max.sample.size",
  "precipitation_daily.cum_inflection2_yearly.mean", "precipitation_daily.cum_inflection2_yearly.sd",
  "precipitation_daily.cum_rate2_yearly.mean",  "precipitation_daily.cum_rate2_yearly.sd",                      "precipitation_daily.cum_n.sample"
)]

simExample <- simExample[, c(
  "X.step.", "T", "RAIN", "solarRadiation", "ET_0", 
  "mean..WATp..of.patches", "mean..ARID..of.patches",
  "mean..biomass..of.patches.with..position.crop.typesOfCrops...0.",
  "mean..biomass..of.patches.with..position.crop.typesOfCrops...1.",
  "mean..yield..of.patches.with..position.crop.typesOfCrops...0.",
  "mean..yield..of.patches.with..position.crop.typesOfCrops...1."
)]

names(simExample) <- c(
  "day", "T", "RAIN", "solarRadiation", "ETr",
  "WATp", "ARID", 
  "wheat_biomass", "rice_biomass", "wheat_yield", "rice_yield"
  )

# create vector with harvest yields per crop and year

years <- 1:floor(nrow(simExample) / 365) - 1

cropYields <- data.frame(
  wheatHarvestDays = (years * 365 + wheatHarvestDay),
  riceHarvestDays = (years * 365 + riceHarvestDay),
  wheat = unique(simExample$wheat_yield),
  rice = unique(simExample$rice_yield)[-1])
cropYields$wheat[1] <- NA # ignore first zero yield in wheat

Show table with parameter values:

parvalues <- c()
for (i in 1:length(simExample_pars))
{
    parvalues <- c(parvalues, simExample_pars[[i]])
}
parvalues <- cbind(names(simExample_pars), parvalues)
knitr::kable(parvalues, 
             format = "html",
             col.names = c("parameter", "values"),
             align = c("l", "r"))
parameter values
solar_annual.min 3
solar_annual.max 7
solar_daily.mean.fluctuation 1
temperature_annual.min.at.2m 15
temperature_annual.max.at.2m 40
temperature_daily.mean.fluctuation 5
temperature_daily.lower.deviation 5
temperature_daily.upper.deviation 5
precipitation_yearly.mean 400
precipitation_yearly.sd 130
precipitation_daily.cum_plateau.value_yearly.mean 0.1
precipitation_daily.cum_inflection1_yearly.mean 40
precipitation_daily.cum_inflection1_yearly.sd 20
precipitation_daily.cum_rate1_yearly.mean 0.15
precipitation_daily.cum_rate1_yearly.sd 0.02
precipitation_daily.cum_plateau.value_yearly.sd 0.05
precipitation_daily.cum_max.sample.size 10
precipitation_daily.cum_inflection2_yearly.mean 200
precipitation_daily.cum_inflection2_yearly.sd 20
precipitation_daily.cum_rate2_yearly.mean 0.05
precipitation_daily.cum_rate2_yearly.sd 0.01
precipitation_daily.cum_n.sample 200


Load function for marking the end of each year:

markEndYears <- function(lengthOfData, offset = 1){
  for (i in 1:lengthOfData)
  {
    if (i %% (365 * offset) == 0)
    {
      abline(v = i, lty = 3)
    }
  }
}

Plot all selected variables:

#--------
# To create a png file:
plotName = "Fig6.png"

fontScale = c(2, 2, 1.2, # c(cex.lab, legend, cex.axis)
              3, 2) # side annotation name and "model"

png(plotName, width = 1000, height = 1200)
#---------
# alternatively, to create eps file:
# plotName = "Fig6.eps"
# 
# fontScale = c(1.8, 1.8, 1.2, # c(cex.lab, legend, cex.axis)
#               2.8, 1.8) # side annotation name and "model"
# 
# extrafont::loadfonts(device = "postscript")
# grDevices::postscript(file = plotName,
#                       pointsize = 10,
#                       width = 10,
#                       height = 12,
#                       horizontal = FALSE,
#                       paper = "special",
#                       onefile = FALSE,
#                       family = "sans",
#                       colormodel = "cmyk")
#---------

layout(matrix(c(1:8, rep(9, 4), rep(10, 2), 11, 12), 
              nrow = 8, ncol = 2), 
       heights = c(1, 1, 1, 1, 1, 1, 1.3, 0.26),
       widths = c(1, 0.15))

yLabs <- c(expression(paste("    Solar\nRadiation (", MJ/m^-2, ")")), 
           "Temperature (ºC)", 
           "Precipitation (mm)", 
           "ETr (mm)", "WATp (%)", "ARID index")
leftMargin <- 7


# 1. Solar radiation
par(mar = c(0.1, leftMargin, 1.1, 0.2), 
    cex.lab = fontScale[1],
    cex.axis = fontScale[3])

plot(1:nrow(simExample), simExample$solarRadiation, type = "l",
     xlab = "", xaxt = 'n',
     ylab = yLabs[1])
markEndYears(nrow(simExample))

# 2. Temperature
plot(1:nrow(simExample), simExample$"T", type = "l",
     xlab = "", xaxt = 'n',
     ylab = yLabs[2])
#lines(1:nrow(simExample), simExample$)
markEndYears(nrow(simExample))

# 3. Precipitation
par(mar = c(1.1, leftMargin, 2.1, 0.2))

barplot(simExample$RAIN, 
        ylab = yLabs[3])
markEndYears(nrow(simExample), offset = 1.2) ; abline(v = 2190 * 1.2, lty = 3)
# not sure why, but barplot() x coordinates do not behave as in plot()

# 4. Reference evapotranspiration
par(mar = c(0.1, leftMargin, 1.1, 0.2))

barplot(simExample$ETr, 
        ylab = yLabs[4])
markEndYears(nrow(simExample), offset = 1.2) ; abline(v = 2190 * 1.2, lty = 3)
# not sure why, but barplot() x coordinates do not behave as in plot()

# 5. Soil water content (%)
plot(1:nrow(simExample), simExample$WATp * 100, 
     xlab = "", xaxt = 'n',
     ylab = yLabs[5],
     type = "l", lwd = 2, ylim = c(0, soil.FC * 110))
abline(h = soil.FC * 100, lty = 2)
abline(h = soil.WP * 100, lty = 2)
markEndYears(nrow(simExample))
text(x = 2000, y = 17.5, labels = "field capacity", 
     cex = fontScale[1], font = 3)
text(x = 2000, y = 4.5, labels = "wilting point", 
     cex = fontScale[1], font = 3)

# 6. ARID index
plot(1:nrow(simExample), simExample$ARID, 
     xlab = "", xaxt = 'n',
     ylab = yLabs[6],
     type = "l", lwd = 2, ylim = c(0, 1))
markEndYears(nrow(simExample))

# 7. Crop biomass and yield
par(mar = c(3.5, leftMargin, 1.1, 0.2))

plot(1:nrow(simExample), simExample$wheat_biomass, type = "l",
     xlab = "", col = "darkred", lwd = 2,
     ylab = expression(paste("Biomass (", g/m^-2, ")")))
lines(1:nrow(simExample), simExample$rice_biomass, col = "darkblue", lwd = 2)
lines(cropYields$wheatHarvestDays, cropYields$wheat, col = "darkred", lwd = 2, lty = 3, pch = 4, cex = 3, type = "o")
lines(cropYields$riceHarvestDays, cropYields$rice, col = "darkblue", lwd = 2, lty = 3, pch = 4, cex = 3, type = "o")
markEndYears(nrow(simExample))

text(x = 2050, y = 700, labels = "harvested", 
     cex = fontScale[1], font = 3)

arrows(x0 = 2050, x1 = 1950, y0 = 620, y1 = 450, length = 0.1, lwd = 2)

mtext("day", side = 1, line = 3, cex = fontScale[3])

CurlyBraces(x0 = 0.02, x1 = 364.08, 
            y0 = -100, y1 = -100, 
            position = 1, 
            bracesSize = 100)
mtext("year (365 days)", side = 1, line = 3, 
      cex = fontScale[3], adj = 0.06, font = 3)

# 8. Legend (crops)
par(mar = c(0,0,0,0), cex = fontScale[2])

plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

legend(x = 0.05, y = 1.1, 
       legend = c('Triticum aestivum (Yecora Rojo)', 'Oryza sativa (IR72)'),
       col = c('darkred', 'darkblue'), lty = c(1, 1), lwd = 2,
       horiz = T, bty = "n")

# 9. Weather model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])

plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

CurlyBraces(x0 = 0.01, x1 = 0.01, 
            y0 = -0.02, y1 = 1.01, 
            position = 4, 
            bracesSize = 0.66,
            label = "Weather\n",
            labelSize = fontScale[4])
mtext(text = "model", 
      side = 4, line = 0.33, cex = fontScale[5])

# 10. Soil Water Balance model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])

plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

CurlyBraces(x0 = 0.01, x1 = 0.01, 
            y0 = -0.02, y1 = 1, 
            position = 4, 
            bracesSize = 0.66,
            label = "Soil Water Balance\n",
            labelSize = fontScale[4])
mtext(text = "model", 
      side = 4, line = 0.33, cex = fontScale[5])

# 11. Crop model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])

plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

CurlyBraces(x0 = 0.01, x1 = 0.01, 
            y0 = 0.1, y1 = 1, 
            position = 4, 
            bracesSize = 0.66,
            label = "Crop\n",
            labelSize = fontScale[4])
mtext(text = "model", 
      side = 4, line = 0.33, cex = fontScale[5])

# 12. empty
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

dev.off()
## png 
##   2
knitr::include_graphics(plotName)

Figure 7. Example of simulated terrains

Load packages:

library(grid)
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 3.6.3
library(png)

Defining terrain random number generator seed (must be available inside “terrains” folder:

listOfSeeds <- c(32, 298, 304, 499)

Setting dimensions:

figScale = 2
terrainDim = 450 * figScale
topHeight = 100 * figScale
bottomHeight = 150 * figScale
leftWidth = 100 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 3 + leftWidth + 2 * internalMargin
figHeight = terrainDim * 4 + topHeight + bottomHeight + 3 * internalMargin

Create margin header images:

headers <- c("seed",
             "Elevation", 
             "Soil texture type",
             "Ecological community\n(cover type)")
namesInFile <- c("seed", "elev", "soil", "ecol")

for (i in 1:length(headers))
{
  if (i == 1) { width = leftWidth } else { width = terrainDim }
  
  png(paste0("Fig7_input/Fig7_", namesInFile[i], ".PNG"), 
    width = width, 
    height = topHeight)

  par(mar = rep(0, 4))

  plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
  
  text(0.5, 0.5, cex = figScale * 3.7,
       labels = headers[i])
  
  dev.off()
}

for (seed in listOfSeeds)
{
  png(paste0("Fig7_input/Fig7_seed=", seed, ".PNG"), 
    width = leftWidth, 
    height = terrainDim)

  par(mar = rep(0, 4))

  plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
  
  text(0.5, 0.5, cex = figScale * 5, adj = .5,
       labels = seed)
  
  dev.off()
}

Create legend images:

# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 6)

png("Fig7_input/Fig7_legend_elev.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = c(2,2,1,2), cex = figScale * 2)

barplot(rep(100, 6), 
        col = rgb(elev - (100/255), elev, 0), 
        axes = FALSE)
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)

dev.off()

soilTextureTypes <- c("Sandy loam", "Loam", 
                      "Sandy clay loam", "Clay loam", 
                      "Sandy clay", "Clay")
soilTextureTypes_colors <- c("#9d6e48", "#eded31", 
                             "#a71b6a", "#345da9",
                             "#7c50a4", "#2d8dbe")
xy <- list(c(-0.1, 1),   c(0.55,  1), 
           c(-0.1, 0.7), c(0.55, 0.7),  
           c(-0.1, 0.4), c(0.55, 0.4))

png("Fig7_input/Fig7_legend_soil.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = rep(0, 4))

plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)

for (i in 1:length(soilTextureTypes))
{
  legend(xy[[i]][1], xy[[i]][2], 
         legend = soilTextureTypes[i], 
         fill = soilTextureTypes_colors[i], 
         bty = "n",  adj = 0.1,
         cex = figScale * 2.5, pt.cex = figScale * 5)
}

dev.off()

ecologicalCommunityTypes <- c("grassland", "wood-grass",
                              "shrubland", "woodland")
ecologicalCommunityTypes_colors <- c("#f16a15", "#eded31",
                                     "#59b03c", "#1d9f78")
xy <- list(c(0, 0.95), c(0.45, 0.95), 
           c(0, 0.5),  c(0.45, 0.5))

png("Fig7_input/Fig7_legend_ecol.PNG", 
    width = terrainDim, 
    height = bottomHeight)

par(mar = rep(0, 4))

plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)

for (i in 1:length(ecologicalCommunityTypes))
{
  legend(xy[[i]][1], xy[[i]][2], 
         legend = ecologicalCommunityTypes[i], 
         fill = ecologicalCommunityTypes_colors[i], 
         bty = "n", adj = 0.1,
         cex = figScale * 2.5, pt.cex = figScale * 5)
}

dev.off()

Create blank image:

png("Fig7_input/Fig7_blank.PNG", 
    width = 10, 
    height = 10)
## Warning in png("Fig7_input/Fig7_blank.PNG", width = 10, height = 10): 'width=10,
## height=10' are unlikely values in pixels
dev.off()
## png 
##   2

Load input images:

elevTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                    listOfSeeds ,".PNG")
soilTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                    listOfSeeds ,"_soilTextureTypes.PNG")
ecolTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", 
                    listOfSeeds ,"_coverTypes.PNG")

# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("Fig7_input/Fig7_seed.PNG",
                 "Fig7_input/Fig7_elev.PNG",
                 "Fig7_input/Fig7_blank.PNG",
                 "Fig7_input/Fig7_soil.PNG",
                 "Fig7_input/Fig7_blank.PNG",
                 "Fig7_input/Fig7_ecol.PNG",
                 # ---
                 paste0("Fig7_input/Fig7_seed=", listOfSeeds[1], ".PNG"),
                 elevTerrains[1],
                 "Fig7_input/Fig7_blank.PNG",
                 soilTerrains[1],
                 "Fig7_input/Fig7_blank.PNG",
                 ecolTerrains[1],
                 # ---
                 rep("Fig7_input/Fig7_blank.PNG", 6),
                 # ---
                 paste0("Fig7_input/Fig7_seed=", listOfSeeds[2], ".PNG"),
                 elevTerrains[2],
                 "Fig7_input/Fig7_blank.PNG",
                 soilTerrains[2],
                 "Fig7_input/Fig7_blank.PNG",
                 ecolTerrains[2],
                 # ---
                 rep("Fig7_input/Fig7_blank.PNG", 6),
                 # ---
                 paste0("Fig7_input/Fig7_seed=", listOfSeeds[3], ".PNG"),
                 elevTerrains[3],
                 "Fig7_input/Fig7_blank.PNG",
                 soilTerrains[3],
                 "Fig7_input/Fig7_blank.PNG",
                 ecolTerrains[3],
                 # ---
                 rep("Fig7_input/Fig7_blank.PNG", 6),
                 # ---
                 paste0("Fig7_input/Fig7_seed=", listOfSeeds[4], ".PNG"),
                 elevTerrains[4],
                 "Fig7_input/Fig7_blank.PNG",
                 soilTerrains[4],
                 "Fig7_input/Fig7_blank.PNG",
                 ecolTerrains[4],
                 # ---
                 "Fig7_input/Fig7_blank.PNG",
                 "Fig7_input/Fig7_legend_elev.PNG",
                 "Fig7_input/Fig7_blank.PNG",
                 "Fig7_input/Fig7_legend_soil.PNG",
                 "Fig7_input/Fig7_blank.PNG",
                 "Fig7_input/Fig7_legend_ecol.PNG"
                 )

grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)

Create Figure 7:

widths = 50 * c(leftWidth/figWidth,
                terrainDim/figWidth,
                internalMargin/figWidth,
                terrainDim/figWidth,
                internalMargin/figWidth,
                terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 50 * c(topHeight/figHeight, 
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        internalMargin/figHeight,
                                        terrainDim/figHeight,
                                        bottomHeight/figHeight) * figScale

lay <- rbind(1:6, 7:12, 13:18, 19:24, 25:30, 31:36, 37:42, 43:48, 49:54)

png("Fig7.png", width = figWidth, height = figHeight)

gridExtra::grid.arrange(grobs = grobs, 
                        layout_matrix = lay,
                        widths = unit(widths, "cm"),
                        heights = unit(heights, "cm"))

dev.off()
## png 
##   2
knitr::include_graphics("Fig7.png")

terrainsData <- read.csv(
    paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfSeeds[1] ,".csv"), skip = 8, nrows = 1)

for (i in 2:length(listOfSeeds))
{
  terrainsData <- rbind(terrainsData, 
                        read.csv(
    paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfSeeds[i] ,".csv"), skip = 8, nrows = 1))
}
  
terrainsData <- as.data.frame(t(terrainsData[, c(89, 18:39, 
                                 94, 96, 99:102, 104:107, 111,
                                 11:16)]))
#names(terrainsData) <- listOfSeeds

Show table with parameter values:

knitr::kable(terrainsData, 
             format = "html",
             col.names = rep("", 4),
             align = c("l", "c", "c", "c"))
randomseed 32 298 304 499
elev_algorithm.style “C#” “C#” “C#” “C#”
elev_featureanglerange 10.355587 13.886082 1.513911 15.282367
elev_inversioniterations 1.32550504 0.61290965 0.96763471 0.05134417
elev_noise 0.4589207 4.4218771 0.9439117 0.7888913
elev_numdepressions 2 9 2 5
elev_numprotuberances 2 1 8 8
elev_numranges 16 52 69 7
elev_numrifts 67 48 61 71
elev_rangeaggregation 0.8162051 0.2617838 0.4092565 0.1425585
elev_rangeheight 36.833480 5.124392 17.353900 30.715513
elev_rangelength 3311 769 1418 2718
elev_riftaggregation 0.59655344 0.04929822 0.18191958 0.73627220
elev_riftheight -30.45546 -22.91505 -43.80874 -32.58051
elev_riftlength 3217 3021 3085 2546
elev_smoothingradius 3.464823 3.464823 3.464823 3.464823
elev_smoothstep 1 1 1 1
elev_valleyaxisinclination 0.7037125 0.7573202 0.0565110 0.6901128
elev_valleyslope 0.011788025 0.002599597 0.004806264 0.010067643
elev_xslope 0.004417135 0.005646562 0.007325109 0.007044657
elev_yslope 0.005514878 0.006371659 0.006696862 0.009834952
flow_do.fill.sinks true true true true
flow_riveraccumulationatstart 34790950 18916696 20365916 43365595
soil_depthnoise 4.645497 40.531347 38.775613 38.513392
soil_formativeerosionrate 1.685376 2.807310 1.731212 1.152460
soil_max.clay 75.65757 61.35114 64.55615 90.23157
soil_max.sand 94.74552 99.87023 97.47498 96.08472
soil_max.silt 99.96852 49.18063 64.34818 95.20315
soil_maxdepth 514.8621 230.2860 387.0962 297.2467
soil_min.clay 16.24237 61.34954 38.91594 73.28332
soil_min.sand 90.25214 99.54846 54.11990 88.72363
soil_min.silt 45.20185 31.31848 59.64769 91.84854
soil_mindepth 241.5896 194.5769 149.2301 147.3834
soil_texturenoise 1.606241 6.412468 1.963127 4.493502
ecol_brushfrequencyinflection 20.10049 14.47833 12.84533 26.76113
ecol_brushfrequencyrate 0.01765001 0.09451887 0.03295388 0.09052492
ecol_grassfrequencyinflection 4.0538838 3.2674253 1.9392423 0.5540513
ecol_grassfrequencyrate 0.008029435 0.133515886 0.047298747 0.109743989
ecol_woodfrequencyinflection 41.68921 54.95257 30.47011 32.99336
ecol_woodfrequencyrate 0.005506242 0.008327215 0.071837624 0.066332965



Video. Example of one-year simulation of a preliminary version of the Integrated Land Unit model (Weather + Soil Water Balance + Land)